class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 4: Procesamiento masivo de datos geoespaciales ### Breve introducción a las Series de tiempo MatÃas Olea <br> <a href="https://orcid.org/0000-0003-0194-7784"> ORCID </a><br> matias.olea@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre 2022</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## ¿Qué aprenderemos en esta unidad? Aprenderemos ... -- 1) Que es una serie de tiempo. -- 2) Extracción de series de tiempo (vector) desde raster temporales. -- 3) Rellenado de datos. -- 4) Remuestreo temporal. -- 5) Componentes de una serie y descomposición. -- 6) Análisis de cambios en tendencias (ciclos). -- 7) AnomalÃas respecto a desviaciones de la media (Z-score). --- ### Series de tiempo Una **serie de tiempo** corresponde a un conjunto de observaciones de variables cuantitativas en el tiempo. Estas variables pueden o no tener distintos comportamiento en el tiempo.
--- ### Series de tiempo Las series las podemos encontrar como estacionarias y no estacionarias .pull-left[ **Estacionaria**. Es aquella cuya media y varianza se mantienen en el tiempo, es decir, son constantes, por lo que facilita más su modelamiento.
] .pull-left[ **No Estacionaria**. Es aquella cuya media y varianza cambian en el tiempo, es decir, no son constantes, por lo que dificulta más su modelamiento.
] --- ### Series de tiempo Las series temporales poseen 4 componentes (3 para algunos autores) los cuales son en simples palabras: -- **Tendencia** (T - trend): corresponde al comportamiento de la media y la varianza en el tiempo. -- **Ciclos** (C - cycle): oscilaciones respecto a la tendencia en largas o medias escalas temporales. Como se desprende de la tendencia, algunos autores no la consideran como un componente propiamente tal. -- **Estacionalidad** (S - seasonality): comportamientos determinados en periodos fijos. Corresponden a cambios dentro de unidades de tiempo de manera regular (como cambios intra-anuales). -- **Aleatorio** (R - random): comportamiento sin un patrón determinado. Es un comportamiento practicamente impredecible. Lo podemos encontrar en la literatura también como "erratico". --- ### Series de tiempo #### Lectura de datos ```r library(terra) library(magrittr) ``` ```r setwd("C:/TuDirectorio/") # Leemos fechas desde 2001 y /10000 para pasar de numeros enteros a decimales mod.ndvi <- rast("MOD13Q1_NDVI_ts492.tif")[[21:492]]/10000 t.dates <- read.csv("MOD13Q1_dates_ts492.csv") # leemos nuestra tabla con fechas dates <- t.dates$dates[21:492] %>% as.Date() # separamos fechas y pasamos a class Date ``` ```r nlyr(mod.ndvi) # número de bandas (1 banda = 1 fecha) ``` ``` ## [1] 472 ``` ```r length(dates) # número de fechas ``` ``` ## [1] 472 ``` --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel ```r xy1 <- cbind(225529, 6064306) %>% as.matrix() # Coord. Plantación con incendio 2017 xy2 <- cbind(230877, 6038355) %>% as.matrix() # Coord. Sucesión forestal xy3 <- cbind(246458, 6078281) %>% as.matrix() # Coord. Bosque nativo deciduo ``` ```r e1 <- extract(mod.ndvi, xy1) %>% as.numeric() e2 <- extract(mod.ndvi, xy2) %>% as.numeric() e3 <- extract(mod.ndvi, xy3) %>% as.numeric() ``` .footnote[ <span style="font-size:9pt"> Como vimos en clases anteriores, en lugar de usar una matriz con las coordenadas XY de un pixel, podemos utilizar un polÃgono incorporando el estadÃstico que queremos usar y si queremos incluir los pixeles que toncan los bordes del polÃno con el argumento touches= </span> ] --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-10-1.png" width="100%" /> --- ### Series de tiempo #### Rellenado de la serie Existen muchos paquetes que nos permiten realizar rellenos o suavizado de nuestros datos como **gapfill**, **oce**, **imputeTS**, **TSstudio**, entre otros. Los más frecuentes para análisis de datos de IV proveniente de datos satelitales es **phenopix** o el software **TIMESAT**. Recientemente se publicó un nuevo paquete llamado **phenofit**<sup>1</sup> que incluye mejoras y optimizaciones a los tÃpicos algoritmos paramétricos para estudiar series de tiempo de vegetación. Para este ejercicio ocuparemos este último paquete. ```r install.packages("phenofit") ``` <center><img src="data:image/png;base64,#phenofit.png" width="800px"/></center> .footnote[ <span style="font-size:9pt"> <sup>1</sup>Kong, D., McVicar, T. R., Xiao, M., Zhang, Y., Peña-Arancibia, J. L., Filippa, G., Xie, Y., Gu, X. (2022). phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods in Ecology and Evolution, 13, 1508– 1527. https://doi.org/10.1111/2041-210X.13870 </span> ] --- ### Series de tiempo #### Rellenado de la serie <center><img src="data:image/png;base64,#phenofit2.png" width="1000px"/></center> .footnote[ <span style="font-size:10pt"> **wHANTS**: El algoritmo HANTS fue desarrollado con base en la transformada de Fourier. Las series de Fourier son sucesiones superpuestas, en un intervalo de tiempo, de una constante con senos y cosenos de creciente múltiplos enteros de la frecuencia original basados en el intervalo de tiempo.</span> <span style="font-size:10pt"> **wSG**: El filtro de Savitzky–Golay corresponde al ajuste de una regresión polinomial en ventanas temporales que se ajustan a los datos de entrada pero "suavizados".</span> <span style="font-size:10pt"> **wWHIT**: es un modelo mas robusto para datos outliers que los otros métodos, y puede captural componentes estacionales de vegetación cuando las series de tiempo tienen más del 30% de los datos de tu IV están contaminados. </span> ] --- ### Series de tiempo #### Rellenado de la serie con Savitzky–Golay ```r library(phenofit) s.sg <- smooth_wSG(y=e1, nptperyear = 23, iters = 3,frame=7) f1 <- s.sg$zs$ziter1 # Resultado de la primera iteración plot(x=dates, y=e1, type="l", ylab="NDVI", xlab="", ylim=c(0,1),lwd=2) lines(x=dates, y=f1, type="l",col="red") ``` <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-13-1.png" width="100%" /> --- ### Series de tiempo #### Remuestreo temporal El remuestreo temporal consiste en cambiar la resolución temporal nativa de nuestros datos, a una menor. Si bien, metodológicamente se pueden realizar downscaling del tiempo, los resultados no mejorarán. --- class: center,middle background-image: url(data:image/png;base64,#tiempo.png) background-size: 90% background-color: black --- ### Series de tiempo #### Remuestreo temporal Para este ejemplo, utilizaremos nuestra serie de tiempo NDVI de MODIS-Terra de 16-dÃas, para remuestrearla temporalmene en datos mensuales. Por lo que de un total de 23 datos al año pasaremos a tener 12. Primero creamos un vector de serie de tiempo (xts) con el paquete **xts** ```r library(xts) xts.e1 <- xts(e1,dates) # serie bruta xts.f1 <- xts(f1,dates) # serie rellenada ``` Con la función **apply.monthly()** remuestreamos nuestros datos a escala mensual con el estadÃstico que escojamos. ```r m.e1 <- apply.monthly(xts.e1,FUN=mean, na.rm=T) m.f1 <- apply.monthly(xts.f1,FUN=mean, na.rm=T) ``` --- ### Series de tiempo #### Remuestreo temporal ```r df_ts <- data.frame(Date=index(m.e1), Raw_monthly=as.numeric(m.e1), Filled_monthly=as.numeric(m.f1)) plot(x=df_ts$Date, y=df_ts$Raw_monthly, type="l", ylab="NDVI", xlab="", main="Monthly Data", ylim=c(0,1), lwd=2) lines(x=df_ts$Date, y=df_ts$Filled_monthly, type="l",col="red") legend("bottomleft", legend=c("Raw monthly", "Filled monthly"), col=c("black", "red"), lty=1, cex=0.8) ``` --- ### Series de tiempo #### Remuestreo temporal <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-17-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ```r library("tseries") ``` **Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test** **H<sub>0</sub>**: La serie de tiempo tiene una tendencia estacionaria. **H<sub>A</sub>**: La serie de tiempo tiene una tendencia no estacionaria. --- ### Series de tiempo #### Test de estacionariedad ```r kpss.test(na.omit(e1), null="Trend",lshort=F) # Serie plantación con Incendio ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: na.omit(e1) ## KPSS Trend = 0.26172, Truncation lag parameter = 17, p-value = 0.01 ``` Se rechaza H<sub>0</sub> por lo tanto la serie es no estacionaria. <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-20-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ```r kpss.test(na.omit(e2), null="Trend",lshort=F) # Serie sucesión forestal ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: na.omit(e2) ## KPSS Trend = 0.24633, Truncation lag parameter = 17, p-value = 0.01 ``` Se rechaza H<sub>0</sub> por lo tanto la serie es no estacionaria. <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-22-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ```r kpss.test(na.omit(e3), null="Trend",lshort=F) # Serie plantación con Incendio ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: na.omit(e3) ## KPSS Trend = 0.047889, Truncation lag parameter = 17, p-value = 0.1 ``` No se rechaza H<sub>0</sub> por lo tanto la serie es estacionaria. <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-24-1.png" width="100%" /> --- ### Series de tiempo #### Descomposición de la serie ```r # TS de la serie mensual bruta ts.rm <- ts(data= df_ts$Raw_monthly, start=c(2001,1), end=c(2021,6), frequency= 12) # TS de la serie mensual rellenada ts.fm <- ts(data= df_ts$Filled_monthly, start=c(2001,1), end=c(2021,6), frequency= 12) ``` ```r dec.raw <- decompose(ts.rm) dec.fil <- decompose(ts.fm) ``` --- ### Series de tiempo #### Descomposición de la serie ```r forecast::autoplot(dec.raw) # Descomposición serie mensual a partir de datos brutos ``` <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-27-1.png" width="100%" /> --- ### Series de tiempo #### Descomposición de la serie ```r forecast::autoplot(dec.fil) # Descomposición serie mensual a partir de datos rellenados ``` <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-28-1.png" width="100%" /> --- ### Series de tiempo #### Análisis de cambios en tendencias (ciclos) ```r library(DBEST) cic <- DBEST(data=ts.rm, data.type = "cyclical", seasonality=12, algorithm = "change detection", breakpoints.no=3, first.level.shift=0.05,second.level.shift=0.1, distance.threshold = "default", duration=24,alpha=0.05) ``` --- ### Series de tiempo #### Análisis de cambios en tendencias (ciclos) ```r plot(cic) ``` <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-30-1.png" width="100%" /> --- ### Series de tiempo #### AnomalÃas respecto a desviaciones de la media (Z-score) ```r df_ts$anom_Zr <- (df_ts$Raw_monthly - mean(df_ts$Raw_monthly)) / sd(df_ts$Raw_monthly) barplot(df_ts$anom_Zr, names=lubridate::year(df_ts$Date), col= ifelse(df_ts$anom_Zr < 0,"red","blue"), ylab="Anomalia") ``` <img src="data:image/png;base64,#Serie-tiempo_files/figure-html/unnamed-chunk-31-1.png" width="100%" /> --- ### BibliografÃa (2012) E. B. Brooks, V. A. Thomas, R. H. Wynne and J. W. Coulston, "Fitting the Multitemporal Curve: A Fourier Series Approach to the Missing Data Problem in Remote Sensing Analysis," in IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 9, pp. 3340-3353. doi: 10.1109/TGRS.2012.2183137. (2004) Jin Chen, Per. Jönsson, Masayuki Tamura, Zhihui Gu, Bunkei Matsushita, Lars Eklundh, A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter, Remote Sensing of Environment, Volume 91, Issues 3–4, (2022) Kong, D., McVicar, T. R., Xiao, M., Zhang, Y., Peña-Arancibia, J. L., Filippa, G., Xie, Y., Gu, X. phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods in Ecology and Evolution, 13, 1508– 1527. https://doi-org.pucv.idm.oclc.org/10.1111/2041-210X.13870 (2015) G. Yang, H. Shen, L. Zhang, Z. He and X. Li, "A Moving Weighted Harmonic Analysis Method for Reconstructing High-Quality SPOT VEGETATION NDVI Time-Series Data," in IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 11, pp. 6008-6021. doi: 10.1109/TGRS.2015.2431315. --- class: inverse middle 